#### ---- Setup ----------------------------------------------------------- ####

suppressMessages({
  library(tidyverse) 
  library(estimatr)
  library(coefplot)
  library(ggthemes)
  library(ggforce)
  library(haven)
  library(lemon)
  library(cowplot)
  library(grid)
  library(cjoint)
})

## Helper functions needed to run the replication code below. 
source("helpers.R")

## Municipal survey of Yonkers residents 
cvs_dat <-
  read_rds("municipal_survey.rds")

## Stacked conjoint dataset (police + residents)
conjoint_dat <-
  read_rds("conjoints_stacked.rds")

## Plotting params
fixed <- 1
spacing <- .1285

## Setup for cjoint (to adjust for restricted randomization in conjoints) 

## Construct the conjoint design for cjoint
attribute_list <- list()
attribute_list[["sex_"]] <- levels(conjoint_dat$sex_)
attribute_list[["age_"]] <- levels(conjoint_dat$age_)
attribute_list[["race_"]] <- levels(conjoint_dat$race_)
attribute_list[["education_"]] <- levels(conjoint_dat$education_)
attribute_list[["residency_"]] <- levels(conjoint_dat$residency_)
attribute_list[["exam_"]] <- levels(conjoint_dat$exam_)
attribute_list[["motivation_"]] <- levels(conjoint_dat$motivation_)
attribute_list[["occupation_"]] <- levels(conjoint_dat$occupation_)

## Specify randomization constraints on education and occupation attributes. 
## Respondents did not see a school teacher or social worker with GED, 
## High School, or Associates education level
constraint_list <- list()
constraint_list[[1]] <- list()
constraint_list[[1]][["education_"]] <- c("GED", "High school", 
                                          "Associates degree")
constraint_list[[1]][["occupation_"]] <- c("School teacher", "Social worker")

## Pass to makeDesign function
cjoint_design <- 
  makeDesign(type = 'constraints',
             attribute.levels = attribute_list,
             constraints = constraint_list)

#### -----Fig. S36: Estimated marginal means for all attributes ----------- ####
community_est_out <- list()
police_est_out <- list()
var_list <- c(
  "sex_",
  "age_",
  "race_",
  "education_",
  "residency_",
  "exam_",
  "motivation_",
  "occupation_"
)

for (i in 1:length(var_list)) {
  col <- sym(var_list[i])
  community_est_out[[i]] <-
    conjoint_dat %>%
    filter(sample == "Civilian sample") %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i],!!col))
  
  police_est_out[[i]] <-
    conjoint_dat %>%
    filter(sample == "Police sample") %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i],!!col))
}

community_est_df <-
  bind_rows(community_est_out) %>%
  ungroup() %>%
  select(Z, estimate, std.error, statistic, p.value, conf.low, conf.high) %>%
  mutate(
    attribute = str_to_title(word(Z, 1, sep = "_")),
    level = word(Z, 2, sep = "_"),
    level = fct_reorder(level, estimate),
    group = "Civilian sample"
  )

police_est_df <-
  bind_rows(police_est_out) %>%
  ungroup() %>%
  select(Z, estimate, std.error, statistic, p.value, conf.low, conf.high) %>%
  mutate(
    attribute = str_to_title(word(Z, 1, sep = "_")),
    level = word(Z, 2, sep = "_"),
    level = fct_reorder(level, estimate),
    group = "Police sample"
  )

gg_df <-
  bind_rows(police_est_df, community_est_df)

plot_dat <-
  gg_df %>%
  mutate(
    attribute = case_when(
      attribute == "Exam" ~ "Performance on civil service exam",
      attribute == "Race" ~ "Race/ethnicity",
      attribute == "Residency" ~ "Length of residency in city",
      attribute == "Occupation" ~ "Previous occupation",
      attribute == "Motivation" ~ "Motivation for becoming a police officer",
      attribute == "Education" ~ "Educational attainment",
      TRUE ~ paste(attribute)
    ),
    attribute = factor(
      attribute,
      levels = c(
        "Race/ethnicity",
        "Sex",
        "Performance on civil service exam",
        "Length of residency in city",
        "Motivation for becoming a police officer",
        "Previous occupation",
        "Educational attainment",
        "Age"
      )
    ),
    group = factor(group, levels = rev(c(
      "Civilian sample",
      "Police sample"
    )))
  )

g <-
  plot_dat %>%
  ggplot(.,
         aes(
           x = estimate,
           y = level,
           group = group,
           color = group,
           fill = group,
           shape = group
         )) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .9),
    size = 1
  ) +
  geom_point(size = 2.5,
             color = "white",
             position = position_dodgev(height = .9)) +
  facet_col(attribute ~ ., scales = "free_y", space = "free") +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_shape_manual("", values = rev(c(22, 21))) +
  labs(x = "Probability of selecting applicant",
       y = "",
       caption = "") +
  theme_tufte(base_size = 14) +
  theme(
    strip.text.x = element_text(hjust = 0, size = 14),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 12),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10,-10,-10,-10)
  )

g

p <-
  g +
  geom_label(
    data = g$data %>%
      filter(level == "Black"),
    aes(label = group),
    hjust = -0.5,
    position = position_dodgev(height = 3.5),
    color = "white",
    size = 3
  )

p

ggsave(p,
       filename = "estimated_mms_all.pdf",
       width = 7,
       height = 10.5)


#### ---- Fig. S37: estimated AMCEs for all attributes -------------------- #### 

## Estimate AMCEs  using all attributes in the design, w/ adjustments for 
## conditionally independent attributes (see Hainmueller at al. 2014 S4.1)
est_community <-
  amce(
    conjoint_chosen ~ sex_ + age_ + race_ + education_ +
      residency_ + exam_ + motivation_ + occupation_,
    design = cjoint_design, 
    cluster = TRUE,
    respondent.id = "person_id",
    data = conjoint_dat %>%
      filter(sample == "Civilian sample" ) %>% 
      select(person_id, conjoint_chosen, sex_, age_, race_, education_,
             residency_, exam_, motivation_, occupation_) %>% 
      na.omit()
  )

est_police <-
  amce(
    conjoint_chosen ~ sex_ + age_ + race_ + education_ +
      residency_ + exam_ + motivation_ + occupation_,
    design = cjoint_design, 
    cluster = TRUE,
    respondent.id = "person_id",
    data = conjoint_dat %>%
      filter(sample == "Police sample") %>% 
      select(person_id, conjoint_chosen, sex_, age_, race_, education_,
             residency_, exam_, motivation_, occupation_) %>% 
      na.omit()
  )

## Combine for plotting:
gg_df <-
  bind_rows(
    summary(est_community)$amce %>% mutate(sample = "Civilian sample" ),
    summary(est_police)$amce %>% mutate(sample = "Police sample")
  ) %>%
  set_names(tolower) %>%
  rename(std.error = `std. err`,
         statistic = `z value`,
         p.value = `pr(>|z|)`) %>%
  mutate(
    sample = factor(sample, levels = rev(c(
      "Civilian sample" , "Police sample"
    )),
    labels = rev(c(
      "Civilian sample", "Police sample"
    ))),
    attribute = str_to_title(word(attribute, 1, sep = "_")),
    conf.low = estimate - 1.96 * std.error,
    conf.high = estimate + 1.96 * std.error,
    attribute_ref = factor(
      attribute,
      levels = c(
        "Race",
        "Sex",
        "Exam",
        "Residency",
        "Motivation",
        "Occupation",
        "Education",
        "Age"
      ),
      labels = c(
        "Race/Ethnicity (reference: White)",
        "Sex (reference: Male)",
        "Civil service exam performance (reference: Top 25%)",
        "Length of residency (reference: Does not live in city)",
        "Motivation for application (reference: Job benefits)",
        "Previous occupation (reference: Police officer in another city)",
        "Education (reference: GED)",
        "Age (reference: 23 years old)"
      )
    ),
    level = factor(
      level,
      levels = c(
        "White",
        "Asian",
        "Hispanic/Latino",
        "Black",
        "Top 25%",
        "Top 15%",
        "Top 10%",
        "Top 5%",
        "Top 1%",
        "Does not live in city",
        "For less than 1 year",
        "For 1-2 years",
        "For 3-5 years",
        "For 6-10 years",
        "For more than 10 years",
        "Job benefits",
        "Job security",
        "Career advancement",
        "Friends/relatives in PD",
        "Excitement of the work",
        "To fight crime",
        "Lifelong dream/aspiration",
        "To help people",
        rev(c("Police officer in another city",
              "Military service",
              "Social worker",
              "Security guard",
              "School teacher",
              "Personal trainer",
              "Construction worker",
              "Server/Bartender",
              "Retail salesperson")),
        "High school",
        "Associates degree",
        "Bachelors degree",
        "Graduate degree",
        "Male",
        "Female",
        "23",
        "25",
        "27",
        "29",
        "31",
        "33",
        "35",
        "37"
      ),
      labels = c(
        "White",
        "Asian",
        "Hispanic/Latino",
        "Black",
        "Top 25%",
        "Top 15%",
        "Top 10%",
        "Top 5%",
        "Top 1%",
        "Does not live in city",
        "Less than 1 year",
        "1-2 years",
        "3-5 years",
        "6-10 years",
        "More than 10 years",
        "Job benefits",
        "Job security",
        "Career advancement",
        "Friends/relatives in PD",
        "Excitement of the work",
        "To fight crime",
        "Lifelong dream/aspiration",
        "To help people",
        rev(c("Police officer in another city",
              "Military service",
              "Social worker",
              "Security guard",
              "School teacher",
              "Personal trainer",
              "Construction worker",
              "Server/Bartender",
              "Retail salesperson")),
        "High school",
        "Associates degree",
        "Bachelors degree",
        "Graduate degree",
        "Male",
        "Female",
        "23",
        "25",
        "27",
        "29",
        "31",
        "33",
        "35",
        "37"
      )
    )
  )

### Plot results 
g <-
  gg_df %>%
  ggplot(.,
         aes(
           x = estimate,
           y = level,
           group = sample,
           color = sample,
           fill = sample,
           shape = sample
         )) +
  geom_vline(
    xintercept = 0,
    color = "black",
    linetype = "dashed",
    size = 0.4
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.9),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.9),
             size = 2.5,
             color = "white") +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_shape_manual("", values = rev(c(22, 21))) +
  labs(x = "Average marginal component effect",
       y = "",
       caption = "") +
  theme_tufte(base_size = 12) +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.length.x = unit(0.2, "cm"),
    axis.ticks.x = element_line(size = 0.4),
    strip.text = element_text(size = 12, hjust = 0),
    axis.line.x = element_line(size = 0.4),
    axis.title.x = element_text(size = 11),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0,
      r = 0.5,
      b = -1,
      l = -0.5
    ), "lines"),
    legend.position = "none"
  )

g

p <-
  g +
  theme(legend.position = "none") + 
  geom_label(
    data = g$data %>%
      filter(level == "Black"),
    aes(label = sample),
    hjust = -.4,
    position = position_dodgev(height = 2.50),
    color = "white",
    size = 3
  )

p

ggsave(p, filename = "estimated_amces_all.pdf", width = 6.5, height = 9)


#### ---- Fig. S38: Estimated marginal means for ordinal outcomes --------- ####
community_est_out <- list()
police_est_out <- list()
var_list <- c(
  "sex_",
  "age_",
  "race_",
  "education_",
  "residency_",
  "exam_",
  "motivation_",
  "occupation_"
)

for (i in 1:length(var_list)) {
  col <- sym(var_list[i])
  community_est_out[[i]] <-
    conjoint_dat %>%
    filter(sample == "Civilian sample" ) %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_rating ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i], !!col))
  
  police_est_out[[i]] <-
    conjoint_dat %>%
    filter(sample == "Police sample") %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_rating ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i], !!col))
}

community_est_df <-
  bind_rows(community_est_out) %>%
  ungroup() %>%
  select(Z, estimate, std.error, statistic, p.value, conf.low, conf.high) %>%
  mutate(
    attribute = str_to_title(word(Z, 1, sep = "_")),
    level = word(Z, 2, sep = "_"),
    level = fct_reorder(level, estimate),
    group = "Civilian sample"
  )

police_est_df <-
  bind_rows(police_est_out) %>%
  ungroup() %>%
  select(Z, estimate, std.error, statistic, p.value, conf.low, conf.high) %>%
  mutate(
    attribute = str_to_title(word(Z, 1, sep = "_")),
    level = word(Z, 2, sep = "_"),
    level = fct_reorder(level, estimate),
    group = "Police sample"
  )

gg_df <-
  bind_rows(police_est_df, community_est_df)

plot_dat <-
  gg_df %>%
  mutate(
    attribute = case_when(
      attribute == "Exam" ~ "Performance on civil service exam",
      attribute == "Race" ~ "Race/ethnicity",
      attribute == "Residency" ~ "Length of residency in city",
      attribute == "Occupation" ~ "Previous occupation",
      attribute == "Motivation" ~ "Motivation for becoming a police officer",
      attribute == "Education" ~ "Educational attainment",
      TRUE ~ paste(attribute)
    ),
    attribute = factor(
      attribute,
      levels = c(
        "Race/ethnicity",
        "Sex",
        "Performance on civil service exam",
        "Length of residency in city",
        "Motivation for becoming a police officer",
        "Previous occupation",
        "Educational attainment",
        "Age"
      )
    ),
    group = factor(group, levels = rev(
      c("Civilian sample",
        "Police sample")
    ))
  )

g <- 
  plot_dat %>%
  ggplot(.,
         aes(
           x = estimate,
           y = level,
           group = group,
           color = group,
           fill = group,
           shape = group
         )) +
  geom_vline(
    xintercept = 4,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .9),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    color = "white",
    position = position_dodgev(height = .9)
  ) +
  facet_col(attribute ~ ., scales = "free_y", space = "free") +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_shape_manual("", values = rev(c(22, 21))) +
  scale_x_continuous(limits = c(4, 6.5)) + 
  labs(x = "Marginal mean on applicant rating scale",
       y = "",
       caption = "") +
  theme_tufte(base_size = 14) +
  theme(
    strip.text.x = element_text(hjust = 0, size = 14),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 12),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 12),
    axis.text.y = element_text(size = 12),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.5,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g

p <-
  g +
  geom_label(
    data = g$data %>%
      filter(level == "Black"
      ),
    aes(label = group),
    hjust = -0.5,
    position = position_dodgev(height = 3.5),
    color = "white",
    size = 3
  )

p

ggsave(p, filename = "estimated_mms_all_ordinal.pdf", 
       width = 7, height = 10.5)

#### ---- Fig. S39: Estimated AMCEs for ordinal outcomes ------------------ ####

## Estimate AMCEs  using all attributes in the design, w/ adjustments for 
## conditionally independent attributes (see Hainmueller at al. 2014 S4.1)
est_community <-
  amce(
    conjoint_rating ~ sex_ + age_ + race_ + education_ +
      residency_ + exam_ + motivation_ + occupation_,
    design = cjoint_design, 
    cluster = TRUE,
    respondent.id = "person_id",
    data = conjoint_dat %>%
      filter(sample == "Civilian sample" ) %>% 
      select(person_id, conjoint_rating, sex_, age_, race_, education_,
             residency_, exam_, motivation_, occupation_) %>% 
      na.omit()
  )

est_police <-
  amce(
    conjoint_rating ~ sex_ + age_ + race_ + education_ +
      residency_ + exam_ + motivation_ + occupation_,
    design = cjoint_design, 
    cluster = TRUE,
    respondent.id = "person_id",
    data = conjoint_dat %>%
      filter(sample == "Police sample") %>% 
      select(person_id, conjoint_rating, sex_, age_, race_, education_,
             residency_, exam_, motivation_, occupation_) %>% 
      na.omit()
  )

## Combine for plotting:
gg_df <-
  bind_rows(
    summary(est_community)$amce %>% mutate(sample = "Civilian sample" ),
    summary(est_police)$amce %>% mutate(sample = "Police sample")
  ) %>%
  set_names(tolower) %>%
  rename(std.error = `std. err`,
         statistic = `z value`,
         p.value = `pr(>|z|)`) %>%
  mutate(
    sample = factor(sample, levels = rev(c(
      "Civilian sample" , "Police sample"
    )),
    labels = rev(c(
      "Civilian sample", "Police sample"
    ))),
    attribute = str_to_title(word(attribute, 1, sep = "_")),
    conf.low = estimate - 1.96 * std.error,
    conf.high = estimate + 1.96 * std.error,
    attribute_ref = factor(
      attribute,
      levels = c(
        "Race",
        "Sex",
        "Exam",
        "Residency",
        "Motivation",
        "Occupation",
        "Education",
        "Age"
      ),
      labels = c(
        "Race/Ethnicity (reference: White)",
        "Sex (reference: Male)",
        "Civil service exam performance (reference: Top 25%)",
        "Length of residency (reference: Does not live in city)",
        "Motivation for application (reference: Job benefits)",
        "Previous occupation (reference: Police officer in another city)",
        "Education (reference: GED)",
        "Age (reference: 23 years old)"
      )
    ),
    level = factor(
      level,
      levels = c(
        "White",
        "Asian",
        "Hispanic/Latino",
        "Black",
        "Top 25%",
        "Top 15%",
        "Top 10%",
        "Top 5%",
        "Top 1%",
        "Does not live in city",
        "For less than 1 year",
        "For 1-2 years",
        "For 3-5 years",
        "For 6-10 years",
        "For more than 10 years",
        "Job benefits",
        "Job security",
        "Career advancement",
        "Friends/relatives in PD",
        "Excitement of the work",
        "To fight crime",
        "Lifelong dream/aspiration",
        "To help people",
        rev(c("Police officer in another city",
              "Military service",
              "Social worker",
              "Security guard",
              "School teacher",
              "Personal trainer",
              "Construction worker",
              "Server/Bartender",
              "Retail salesperson")),
        "High school",
        "Associates degree",
        "Bachelors degree",
        "Graduate degree",
        "Male",
        "Female",
        "23",
        "25",
        "27",
        "29",
        "31",
        "33",
        "35",
        "37"
      )
    )
  )

g <-
  gg_df %>%
  ggplot(.,
         aes(
           x = estimate,
           y = level,
           group = sample,
           color = sample,
           fill = sample,
           shape = sample
         )) +
  geom_vline(
    xintercept = 0,
    color = "black",
    linetype = "dashed",
    size = 0.4
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = 0.9),
    size = 1
  ) +
  geom_point(position = position_dodgev(height = 0.9),
             size = 2.5,
             color = "white") +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_shape_manual("", values = rev(c(22, 21))) +
  labs(x = "Average marginal component effect",
       y = "",
       caption = "") +
  theme_tufte(base_size = 12) +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.length.x = unit(0.2, "cm"),
    axis.ticks.x = element_line(size = 0.4),
    strip.text = element_text(size = 12, hjust = 0),
    axis.line.x = element_line(size = 0.4),
    axis.title.x = element_text(size = 11),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0,
      r = 0.5,
      b = -1,
      l = -0.5
    ), "lines"),
    legend.position = "none"
  )

g

p <-
  g +
  theme(legend.position = "none") + 
  geom_label(
    data = g$data %>%
      filter(level == "Black"),
    aes(label = sample),
    hjust = -0.8,
    position = position_dodgev(height = 2.5),
    color = "white",
    size = 3
  )

p

ggsave(p, filename = "estimated_amces_all_ordinal.pdf", 
       width = 6.5, height = 9)


#### ---- Fig. S40: Estimated conditional means by race and sex ----------- ####
gg_df <-
  conjoint_dat %>%
  group_by(race_, sex_, sample) %>%
  do(tidy(lm_robust(
    conjoint_chosen ~ 1, clusters = person_id, data = .
  ))) %>%
  ungroup() %>%
  mutate(sex_ = factor(
    sex_,
    levels = c("Male", "Female"),
    labels = c("Male applicants:", "Female applicants:")
  ),
  X = factor(race_)) %>%
  select(X,
         sample,
         estimate,
         std.error,
         statistic,
         p.value,
         conf.low,
         conf.high,
         sex_,
         race_) %>%
  bind_rows(
    ## Differences between Male and Female applicants
    conjoint_dat %>%
      group_by(sample, race_) %>%
      do(tidy(
        lm_robust(conjoint_chosen ~ sex_,
                  clusters = person_id,
                  data = .)
      )) %>%
      ungroup() %>%
      filter(term != "(Intercept)")  %>%
      mutate(X = "Female v. Male") %>%
      select(
        X,
        sample,
        estimate,
        std.error,
        statistic,
        p.value,
        conf.low,
        conf.high,
        race_
      )
  )

## Plot for male applicants (reference category)
g1 <- 
  gg_df %>%
  filter(sex_ == "Male applicants:") %>%
  ggplot(.,
         aes(x = estimate,
             y = X,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(
    expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 

## Plot for female applicants
g2 <- 
  gg_df %>%
  filter(sex_ == "Female applicants:") %>%
  ggplot(.,
         aes(x = estimate,
             y = X,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(
    expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 

## Differences 
g3 <- 
  gg_df %>%
  filter(X == "Female v. Male") %>%
  ggplot(.,
         aes(x = estimate,
             y = race_,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(
    expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 


g <- 
  egg::ggarrange(
    g2 + labs(x = "Conditional probability of\nselection: Female applicant") + 
      theme(strip.text.x = element_text(size = 10, color = "white"),
            strip.placement = "outside",
            panel.border = element_rect(fill = NA, colour = "black")) +
      geom_label(
        data = g2$data %>%
          filter(race_ == "White") %>% 
          mutate(estimate = 0.52),
        aes(label = sample),
        hjust = -.01,
        position = position_dodgev(height = 0.45),
        color = "white",
        size = 2.5
      ),
    g1 + labs(x = "Conditional probability of\nselection: Male applicant") + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference: average effect of\nFemale (v. Male) applicant")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )


## Convert to grob and export
ggsave(
  filename = "conjoint_race_sex_means.pdf",
  ggplotify::as.ggplot(g),
  width = 7.5,
  height = 2 + fixed + spacing * (nrow(g1$data))
)

#### ---- Fig. S41: Estimated conditional means by race and exam ---------- ####
gg_df <- 
  conjoint_dat %>% 
  group_by(race_, exam_, sample) %>% 
  do(tidy(lm_robust(conjoint_chosen ~ 1, clusters = person_id, data = .))) %>% 
  ungroup() %>% 
  mutate(
    race_ = factor(race_, 
                   levels = c("White", "Black", "Hispanic/Latino", "Asian"),
                   labels = c("White applicants:", "Black applicants:", 
                              "Hispanic/Latino applicants:", 
                              "Asian applicants:")),
    X = race_
  ) %>% 
  select(X,
         sample,
         estimate,
         std.error,
         statistic,
         p.value,
         conf.low,
         conf.high,
         race_,
         exam_) %>%
  bind_rows(
    ## Between sample differences
    conjoint_dat %>%
      group_by(exam_, race_) %>%
      do(tidy(
        lm_robust(conjoint_chosen ~ sample, clusters = person_id, data = .)
      )) %>%
      ungroup() %>%
      filter(term != "(Intercept)") %>%
      mutate(
        race_ = factor(race_, 
                       levels = c("White", "Black", "Hispanic/Latino", "Asian"),
                       labels = c("White applicants:", "Black applicants:", 
                                  "Hispanic/Latino applicants:", 
                                  "Asian applicants:")),
        X = "Civilian v. Police") %>%
      select(
        X,
        estimate,
        std.error,
        statistic,
        p.value,
        conf.low,
        conf.high,
        race_,
        exam_
      )
  ) 

## Plot results for civilian sample
g1 <- 
  gg_df %>%
  filter(sample == "Civilian sample") %>%
  ggplot(data = .,
         aes(
           y = exam_,
           x = estimate,
           fill = exam_,
           color = exam_,
           group = exam_
         )) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.25
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    size = 1,
    height = 0
  ) +
  geom_point(
    pch = 22,
    color = "white",
    size = 2.5
  ) +
  facet_col( ~ race_) +
  scale_color_grey(start = 0.2, end = 0.8) +
  scale_fill_grey(start = 0.2, end = 0.8) +
  labs(y = "",
       x = "Probability of selection",
       caption = "") +
  theme_tufte(base_size = 12) +
  theme(
    strip.text.x = element_text(hjust = 0, size = 12),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.2, color = "black"),
    axis.title.x = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = -0.5
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g1

## Repeat for police sample
g2 <- 
  gg_df %>%
  filter(sample == "Police sample") %>%
  ggplot(data = .,
         aes(
           y = exam_,
           x = estimate,
           fill = exam_,
           color = exam_,
           group = exam_
         )) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.25
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    size = 1,
    height = 0
  ) +
  geom_point(
    pch = 21,
    color = "white",
    size = 2.5
  ) +
  facet_col( ~ race_) +
  scale_color_grey(start = 0.2, end = 0.8) +
  scale_fill_grey(start = 0.2, end = 0.8) +
  labs(y = "",
       x = "Probability of selection",
       caption = "") +
  theme_tufte(base_size = 12) +
  theme(
    strip.text.x = element_text(hjust = 0, size = 12),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.2, color = "black"),
    axis.title.x = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g2

## Between sample differences
g3 <- 
  gg_df %>%
  filter(X == "Civilian v. Police") %>%
  ggplot(data = .,
         aes(
           y = exam_,
           x = estimate,
           fill = exam_,
           color = exam_,
           group = exam_
         )) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.25
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    size = 1,
    height = 0
  ) +
  geom_point(
    pch = 21,
    color = "white",
    size = 2.5
  ) +
  facet_col( ~ race_) +
  scale_color_grey(start = 0.2, end = 0.8) +
  scale_fill_grey(start = 0.2, end = 0.8) +
  scale_x_continuous(
    expand = c(0.01, 0.01)) +
  labs(y = "",
       x = "Difference between Civilian and Police sample",
       caption = "") +
  theme_tufte(base_size = 12) +
  theme(
    strip.text.x = element_text(hjust = 0, size = 12),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.2, color = "black"),
    axis.title.x = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g3


p <-
  g1 +
  theme(legend.position = "none") + 
  geom_label(
    data = g1$data %>%
      filter(race_ == "White applicants:"),
    aes(label = exam_),
    nudge_x = 0.25,
    color = "white",
    size = 2.5
  )

p


g <- 
  egg::ggarrange(
    p + labs(x = "Conditional probability of\nselection: civilian sample") + 
      scale_x_continuous(limits = c(0.15, 0.80), 
                         expand = c(0.01, 0.01)) +
      theme(panel.border = element_rect(fill = NA, colour = "black")),
    g2 + labs(x = "Conditional probability of\nselection: police sample") +
      scale_x_continuous(limits = c(0.15, 0.80), 
                         expand = c(0.01, 0.01)) +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Between sample difference:\ncivilian v. police sample") + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )

ggsave(g, filename = "conjoint_race_exam_means.pdf", 
       width = 7.5, height = 8)

#### ---- Differences in conditional means for race/ethnicity by exam ----- #### 

## Conditional marginal means for white applicants
cmean_white_est <- 
  conjoint_dat %>%
  filter(race_ == "White") %>%
  group_by(sample, exam_) %>%
  do(tidy(lm_robust(
    conjoint_chosen ~ 1,
    clusters = person_id,
    data = .
  ))) %>%
  mutate(X = "White")

## Conditional marginal means for black applicants
cmean_black_est <-
  conjoint_dat %>%
  filter(race_ == "Black") %>%
  group_by(sample, exam_) %>%
  do(tidy(lm_robust(
    conjoint_chosen ~ 1,
    clusters = person_id,
    data = .
  ))) %>%
  mutate(X = "Black")

## Conditional marginal means for hispanic/latino applicants
cmean_latin_est <-
  conjoint_dat %>%
  filter(race_ == "Hispanic/Latino") %>%
  group_by(sample, exam_) %>%
  do(tidy(lm_robust(
    conjoint_chosen ~ 1,
    clusters = person_id,
    data = .
  ))) %>%
  mutate(X = "Hispanic/Latino")

## Conditional marginal means for asian applicants
cmean_asian_est <-
  conjoint_dat %>%
  filter(race_ == "Asian") %>%
  group_by(sample, exam_) %>%
  do(tidy(lm_robust(
    conjoint_chosen ~ 1,
    clusters = person_id,
    data = .
  ))) %>%
  mutate(X = "Asian")


## Differences between black and white applicants
cmean_wb_diff_est <- 
  conjoint_dat %>%
  filter(race_ %in% c("Black", "White")) %>%
  mutate(X = as.numeric(race_ != "White")) %>%
  group_by(sample, exam_) %>%
  do(tidy(
    lm_robust(
      conjoint_chosen ~ X,
      clusters = person_id,
      data = .
    )
  )) %>% 
  filter(term != "(Intercept)")  %>% 
  mutate(X = "Black v. White")

## Differences between hispanic/latino and white applicants
cmean_wh_diff_est <- 
  conjoint_dat %>%
  filter(race_ %in% c("Hispanic/Latino", "White")) %>%
  mutate(X = as.numeric(race_ != "White")) %>%
  group_by(sample, exam_) %>%
  do(tidy(
    lm_robust(
      conjoint_chosen ~ X,
      clusters = person_id,
      data = .
    )
  )) %>% 
  filter(term != "(Intercept)")  %>% 
  mutate(X = "Hispanic/Latino v. White")

## Differences between asian and white applicants
cmean_wa_diff_est <- 
  conjoint_dat %>%
  filter(race_ %in% c("Asian", "White")) %>%
  mutate(X = as.numeric(race_ != "White")) %>%
  group_by(sample, exam_) %>%
  do(tidy(
    lm_robust(
      conjoint_chosen ~ X,
      clusters = person_id,
      data = .
    )
  )) %>% 
  filter(term != "(Intercept)")  %>% 
  mutate(X = "Asian v. White")

gg_cmean_est <-
  bind_rows(
    cmean_white_est,
    cmean_black_est,
    cmean_latin_est,
    cmean_asian_est,
    cmean_wb_diff_est,
    cmean_wh_diff_est,
    cmean_wa_diff_est
  ) %>%
  ungroup() %>%
  select(sample,
         X,
         exam_,
         estimate,
         std.error,
         statistic,
         p.value,
         conf.low,
         conf.high)

## Plot for white applicants (reference category)
g1 <- 
  gg_cmean_est %>%
  filter(X == "White") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 

#### -------- Fig. S42: Export plot for Black v. White -------------------- ####
g2 <- 
  gg_cmean_est %>%
  filter(X == "Black") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 

g3 <- 
  gg_cmean_est %>%
  filter(X == "Black v. White") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_shape_manual("", values = rev(c(22, 21))) +
  scale_x_continuous(expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.5,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g_black <- 
  egg::ggarrange(
    g2 + labs(x = "Conditional probability of\nselection: Black applicant") + 
      theme(strip.text.x = element_text(size = 10, color = "white"),
            strip.placement = "outside",
            panel.border = element_rect(fill = NA, colour = "black")) +
      geom_label(
        data = g2$data %>%
          filter(exam_ == "Top 25%") %>% 
          mutate(estimate = 0.6),
        aes(label = sample),
        hjust = -.01,
        position = position_dodgev(height = 0.5),
        color = "white",
        size = 2.5
      ),
    g1 + labs(x = "Conditional probability of\nselection: White applicant") + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference: average effect of\nBlack (v. White) applicant")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )


## Convert to grob and export
ggsave(
  filename = "conjoint_interactions_black.pdf",
  ggplotify::as.ggplot(g_black),
  width = 7.5, 
  height = 2 + fixed + spacing *(nrow(g1$data))
)


#### -------- Fig. S43: Export plot for Hispanic/Latino v. White ---------- ####
g2 <- 
  gg_cmean_est %>%
  filter(X == "Hispanic/Latino") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 

g3 <- 
  gg_cmean_est %>%
  filter(X == "Hispanic/Latino v. White") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_shape_manual("", values = rev(c(22, 21))) +
  scale_x_continuous(expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.5,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g_latino <- 
  egg::ggarrange(
    g2 + labs(x = "Conditional probability of\nselection: Hispanic/Latino applicant") + 
      theme(strip.text.x = element_text(size = 10, color = "white"),
            strip.placement = "outside",
            panel.border = element_rect(fill = NA, colour = "black")) +
      geom_label(
        data = g2$data %>%
          filter(exam_ == "Top 25%") %>% 
          mutate(estimate = 0.55),
        aes(label = sample),
        hjust = -.01,
        position = position_dodgev(height = 0.5),
        color = "white",
        size = 2.5
      ),
    g1 + labs(x = "Conditional probability of\nselection: White applicant") + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference: average effect of\nHispanic/Latino (v. White) applicant")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )


## Convert to grob and export
ggsave(
  filename = "conjoint_interactions_latino.pdf",
  ggplotify::as.ggplot(g_latino),
  width = 7.5, 
  height = 2 + fixed + spacing *(nrow(g1$data))
)


#### -------- Fig. S44: Export plot for Asian v. White -------------------- ####
g2 <- 
  gg_cmean_est %>%
  filter(X == "Asian") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 

g3 <- 
  gg_cmean_est %>%
  filter(X == "Asian v. White") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_shape_manual("", values = rev(c(22, 21))) +
  scale_x_continuous(expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.5,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g_asian <- 
  egg::ggarrange(
    g2 + labs(x = "Conditional probability of\nselection: Asian applicant") + 
      theme(strip.text.x = element_text(size = 10, color = "white"),
            strip.placement = "outside",
            panel.border = element_rect(fill = NA, colour = "black")) +
      geom_label(
        data = g2$data %>%
          filter(exam_ == "Top 25%") %>% 
          mutate(estimate = 0.55),
        aes(label = sample),
        hjust = -.01,
        position = position_dodgev(height = 0.5),
        color = "white",
        size = 2.5
      ),
    g1 + labs(x = "Conditional probability of\nselection: White applicant") + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference: average effect of\nAsian (v. White) applicant")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )


## Convert to grob and export
ggsave(
  filename = "conjoint_interactions_asian.pdf",
  ggplotify::as.ggplot(g_asian),
  width = 7.5, 
  height = 2 + fixed + spacing *(nrow(g1$data))
)






#### -------- Fig. S45: Differences for white v. non-white pairs ---------- ####


## Let's create identifier for the pairs, and also flag observations that
## involve pairwise comparisons between USA and another country 
npairs <- nrow(conjoint_dat)/2
paired_conjoint <- 
  conjoint_dat %>% 
  arrange(person_id, profile) %>% 
  mutate(pair_id = rep(1:npairs, each = 2)) %>% 
  group_by(pair_id) %>% 
  filter(!is.na(conjoint_chosen)) %>% 
  ungroup()

## First profile
gg1 <-
  left_join(
    paired_conjoint %>% filter(profile == "p_cand_1_1"),
    paired_conjoint %>% filter(profile == "p_cand_1_2"),
    by = c("person_id", "pair_id")
  ) %>%
  filter(race_.x != race_.y) %>% 
  mutate(flag = case_when(
    race_.x == "White" | race_.y == "White" ~ 1, 
    TRUE ~ 0
  )) %>% 
  filter(flag == 1)

## Second profile
gg2 <-
  left_join(
    paired_conjoint %>% filter(profile == "p_cand_2_1"),
    paired_conjoint %>% filter(profile == "p_cand_2_2"),
    by = c("person_id", "pair_id")
  ) %>%
  filter(race_.x != race_.y) %>% 
  mutate(flag = case_when(
    race_.x == "White" | race_.y == "White" ~ 1, 
    TRUE ~ 0
  )) %>% 
  filter(flag == 1)

## Third profile
gg3 <-
  left_join(
    paired_conjoint %>% filter(profile == "p_cand_3_1"),
    paired_conjoint %>% filter(profile == "p_cand_3_2"),
    by = c("person_id", "pair_id")
  ) %>%
  filter(race_.x != race_.y) %>% 
  mutate(flag = case_when(
    race_.x == "White" | race_.y == "White" ~ 1, 
    TRUE ~ 0
  )) %>% 
  filter(flag == 1)

## Fourth profile
gg4 <-
  left_join(
    paired_conjoint %>% filter(profile == "p_cand_4_1"),
    paired_conjoint %>% filter(profile == "p_cand_4_2"),
    by = c("person_id", "pair_id")
  ) %>%
  filter(race_.x != race_.y) %>% 
  mutate(flag = case_when(
    race_.x == "White" | race_.y == "White" ~ 1, 
    TRUE ~ 0
  )) %>% 
  filter(flag == 1)

## Fifth profile
gg5 <-
  left_join(
    paired_conjoint %>% filter(profile == "p_cand_5_1"),
    paired_conjoint %>% filter(profile == "p_cand_5_2"),
    by = c("person_id", "pair_id")
  ) %>%
  filter(race_.x != race_.y) %>% 
  mutate(flag = case_when(
    race_.x == "White" | race_.y == "White" ~ 1, 
    TRUE ~ 0
  )) %>% 
  filter(flag == 1)

these <- c(paste(gg1$person_id, gg1$pair_id),
           paste(gg2$person_id, gg2$pair_id),
           paste(gg3$person_id, gg3$pair_id),
           paste(gg4$person_id, gg4$pair_id),
           paste(gg5$person_id, gg5$pair_id))

any(duplicated(these)) # Sanity check

paired_conjoint <- 
  paired_conjoint %>% 
  mutate(index = paste(person_id, pair_id),
         white_pair = as.numeric(index %in% these)) 


paired_white_est <- 
  paired_conjoint %>% 
  filter(white_pair == 1) %>% 
  filter(race_ == "White") %>%
  group_by(sample, exam_) %>%
  do(tidy(lm_robust(
    conjoint_chosen ~ 1,
    clusters = person_id,
    data = .
  ))) %>%
  mutate(X = "White")


paired_nwhite_est <- 
  paired_conjoint %>% 
  filter(white_pair == 1) %>% 
  filter(race_ != "White") %>%
  group_by(sample, exam_) %>%
  do(tidy(lm_robust(
    conjoint_chosen ~ 1,
    clusters = person_id,
    data = .
  ))) %>%
  mutate(X = "Non-White")

paired_diff_est <- 
  paired_conjoint %>% 
  filter(white_pair == 1) %>% 
  mutate(X = as.numeric(race_ != "White")) %>%
  group_by(sample, exam_) %>%
  do(tidy(lm_robust(
    conjoint_chosen ~ X,
    clusters = person_id,
    data = .
  ))) %>%
  filter(term != "(Intercept)") %>% 
  mutate(X = "Non-White v. White")

gg_paired_est <-
  bind_rows(paired_white_est,
            paired_nwhite_est,
            paired_diff_est) %>%
  ungroup() %>%
  select(sample,
         X,
         exam_,
         estimate,
         std.error,
         statistic,
         p.value,
         conf.low,
         conf.high)

## Plot for white applicants (reference category)
g1 <- 
  gg_paired_est %>%
  filter(X == "White") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 

g2 <- 
  gg_paired_est %>%
  filter(X == "Non-White") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 

g3 <- 
  gg_paired_est %>%
  filter(X == "Non-White v. White") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 

g_paired <- 
  egg::ggarrange(
    g2 + labs(x = "Conditional probability of\nselection: Non-White applicant") + 
      theme(strip.text.x = element_text(size = 10, color = "white"),
            strip.placement = "outside",
            panel.border = element_rect(fill = NA, colour = "black")) +
      geom_label(
        data = g2$data %>%
          filter(exam_ == "Top 25%") %>% 
          mutate(estimate = 0.6),
        aes(label = sample),
        hjust = -.01,
        position = position_dodgev(height = 0.5),
        color = "white",
        size = 2.5
      ),
    g1 + labs(x = "Conditional probability of\nselection: White applicant") + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference: average effect of\nNon-White (v. White) applicant")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )


## Convert to grob and export
ggsave(
  filename = "conjoint_interactions_paired.pdf",
  ggplotify::as.ggplot(g_paired),
  width = 7.5, 
  height = 2 + fixed + spacing *(nrow(g1$data))
)

#### ---- Fig. S46: Estimated conditional means by exam and sex ----------- ####
gg_df <-
  conjoint_dat %>%
  group_by(exam_, sex_, sample) %>%
  do(tidy(lm_robust(
    conjoint_chosen ~ 1, clusters = person_id, data = .
  ))) %>%
  ungroup() %>%
  mutate(sex_ = factor(
    sex_,
    levels = c("Male", "Female"),
    labels = c("Male applicants:", "Female applicants:")
  ),
  X = factor(exam_)) %>%
  select(X,
         sample,
         estimate,
         std.error,
         statistic,
         p.value,
         conf.low,
         conf.high,
         sex_,
         exam_) %>%
  bind_rows(
    ## Between sample differences
    conjoint_dat %>%
      group_by(exam_, sex_) %>%
      do(tidy(
        lm_robust(conjoint_chosen ~ sample, clusters = person_id, data = .)
      )) %>%
      ungroup() %>%
      filter(term != "(Intercept)") %>%
      mutate(sex_ = factor(
        sex_,
        levels = c("Male", "Female"),
        labels = c("Male applicants:", "Female applicants:")
      ),
      X = "Civilian v. Police") %>%
      select(
        X,
        estimate,
        std.error,
        statistic,
        p.value,
        conf.low,
        conf.high,
        sex_,
        exam_
      )
  ) %>%
  bind_rows(
    ## Differences between Male and Female applicants
    conjoint_dat %>%
      group_by(sample, exam_) %>%
      do(tidy(
        lm_robust(conjoint_chosen ~ sex_,
                  clusters = person_id,
                  data = .)
      )) %>%
      ungroup() %>%
      filter(term != "(Intercept)")  %>%
      mutate(X = "Female v. Male") %>%
      select(
        X,
        sample,
        estimate,
        std.error,
        statistic,
        p.value,
        conf.low,
        conf.high,
        exam_
      )
  )


## Plot results for civilian sample
g1 <- 
  gg_df %>%
  filter(sample == "Civilian sample", X != "Female v. Male") %>%
  ggplot(data = .,
         aes(
           y = exam_,
           x = estimate,
           fill = exam_,
           color = exam_,
           group = exam_
         )) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.25
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    size = 1,
    height = 0
  ) +
  geom_point(
    pch = 22,
    color = "white",
    size = 2.5
  ) +
  facet_col( ~ sex_) +
  scale_color_grey(start = 0.2, end = 0.8) +
  scale_fill_grey(start = 0.2, end = 0.8) +
  scale_x_continuous(limits = c(0.28, 0.75), 
                     expand = c(0.01, 0.01)) +
  labs(y = "",
       x = "Probability of selection",
       caption = "") +
  theme_tufte(base_size = 12) +
  theme(
    strip.text.x = element_text(hjust = 0, size = 12),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.2, color = "black"),
    axis.title.x = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = -0.5
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g1


## Repeat for police sample
g2 <- 
  gg_df %>%
  filter(sample == "Police sample", X != "Female v. Male") %>%
  ggplot(data = .,
         aes(
           y = exam_,
           x = estimate,
           fill = exam_,
           color = exam_,
           group = exam_
         )) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.25
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    size = 1,
    height = 0
  ) +
  geom_point(
    pch = 21,
    color = "white",
    size = 2.5
  ) +
  facet_col( ~ sex_) +
  scale_color_grey(start = 0.2, end = 0.8) +
  scale_fill_grey(start = 0.2, end = 0.8) +
  scale_x_continuous(limits = c(0.28, 0.75), 
                     expand = c(0.01, 0.01)) +
  labs(y = "",
       x = "Probability of selection",
       caption = "") +
  theme_tufte(base_size = 12) +
  theme(
    strip.text.x = element_text(hjust = 0, size = 12),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.2, color = "black"),
    axis.title.x = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g2


## Between sample differences
g3 <- 
  gg_df %>%
  filter(X == "Civilian v. Police") %>%
  ggplot(data = .,
         aes(
           y = exam_,
           x = estimate,
           fill = exam_,
           color = exam_,
           group = exam_
         )) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.25
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    size = 1,
    height = 0
  ) +
  geom_point(
    pch = 21,
    color = "white",
    size = 2.5
  ) +
  facet_col( ~ sex_) +
  scale_color_grey(start = 0.2, end = 0.8) +
  scale_fill_grey(start = 0.2, end = 0.8) +
  scale_x_continuous(#limits = c(0.28, 0.75), 
    expand = c(0.01, 0.01)) +
  labs(y = "",
       x = "Difference between Civilian and Police sample",
       caption = "") +
  theme_tufte(base_size = 12) +
  theme(
    strip.text.x = element_text(hjust = 0, size = 12),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.2, color = "black"),
    axis.title.x = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g3


p <-
  g1 +
  theme(legend.position = "none") + 
  geom_label(
    data = g1$data %>%
      filter(sex_ == "Male applicants:"),
    aes(label = exam_),
    nudge_x = 0.1,
    color = "white",
    size = 2.5
  )

p

g <- 
  egg::ggarrange(
    p + labs(x = "Conditional probability of\nselection: civilian sample") + 
      theme(
        panel.border = element_rect(fill = NA, colour = "black")),
    g2 + labs(x = "Conditional probability of\nselection: police sample") +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Between sample difference:\ncivilian v. police sample") + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )


ggsave(g, 
       filename = "conjoint_exam_sex_means.pdf",
       width = 7.5, 
       height = 5)


#### ---- Fig. S47: Differences in conditional means for sex by exam ------ ####

## Plot for male applicants (reference category)
scores <- c("Top 25%", "Top 15%", "Top 10%", "Top 5%", "Top 1%")

g1 <- 
  gg_df %>%
  filter(sex_ == "Male applicants:", X %in% scores) %>%
  mutate(X = factor(X, levels = scores)) %>% 
  ggplot(.,
         aes(x = estimate,
             y = X,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous(
    expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 

## Plot for female applicants
g2 <- 
  gg_df %>%
  filter(sex_ == "Female applicants:", X %in% scores) %>%
  mutate(X = factor(X, levels = scores)) %>% 
  ggplot(.,
         aes(x = estimate,
             y = X,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous( 
    expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 



## Differences 
g3 <- 
  gg_df %>%
  filter(X == "Female v. Male") %>%
  ggplot(.,
         aes(x = estimate,
             y = exam_,
             color = sample,
             group = sample,
             fill = sample,
             shape = sample)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    position = position_dodgev(height = .5)
  ) +
  scale_color_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_fill_manual("", values = rev(c("#000000", "#0072B2"))) +
  scale_x_continuous( 
    expand = c(0.01, 0.01)) +
  labs(x = "",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "none",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 


g <- 
  egg::ggarrange(
    g2 + labs(x = "Conditional probability of\nselection: Female applicant") + 
      theme(strip.text.x = element_text(size = 10, color = "white"),
            strip.placement = "outside",
            panel.border = element_rect(fill = NA, colour = "black")) +
      geom_label(
        data = g2$data %>%
          filter(exam_ == "Top 25%") %>% 
          mutate(estimate = 0.52),
        aes(label = sample),
        hjust = -.01,
        position = position_dodgev(height = 0.45),
        color = "white",
        size = 2.5
      ),
    g1 + labs(x = "Conditional probability of\nselection: Male applicant") + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference: average effect of\nFemale (v. Male) applicant")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )

g

## Convert to grob and export
ggsave(
  filename = "conjoint_interactions_sex.pdf",
  ggplotify::as.ggplot(g),
  width = 7.5, 
  height = 2 + fixed + spacing *(nrow(g1$data))
)


#### ---- Heterogeneity by respondent background characteristics ---------- ####

## Merge civilian subset of conjoint experiment with survey data
conjoint_merged <-
  cvs_dat %>%
  select(person_id, y_ltc_idx_t0, starts_with("x_")) %>%
  right_join(conjoint_dat %>%
               filter(sample == "Civilian sample"),
             by = "person_id")

#### ---- Fig. S48: Differences in marginal means by partisanship --------- ####
community_est_dem0 <- list()
community_est_dem1 <- list()
community_est_dem <- list()

var_list <- c(
  "sex_",
  "race_",
  "residency_",
  "exam_",
  "motivation_"
)

for (i in 1:length(var_list)) {
  col <- sym(var_list[i])
  community_est_dem0[[i]] <-
    conjoint_merged %>%
    filter(x_pid_3_t0 != "Democrat") %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i], !!col),
           X = "No")
  
  community_est_dem1[[i]] <-
    conjoint_merged %>%
    filter( x_pid_3_t0 == "Democrat") %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i], !!col),
           X = "Yes")
  
  community_est_dem[[i]] <-
    conjoint_merged %>%
    mutate(diff = as.numeric(x_pid_3_t0 == "Democrat")) %>% 
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ diff, data = .,
      clusters = person_id
    ))) %>%
    filter(term != "(Intercept)") %>% 
    mutate(Z = paste0(var_list[i], !!col),
           X = "Difference")
  
}


gg_community_mm_pid <-
  bind_rows(community_est_dem0,
            community_est_dem1,
            community_est_dem) %>%
  ungroup() %>%
  select(Z, X, estimate, std.error, statistic, p.value, conf.low, 
         conf.high) %>%
  mutate(
    attribute = str_to_title(word(Z, 1, sep = "_")),
    attribute_ref = factor(
      attribute,
      levels = c("Race",
                 "Sex",
                 "Exam",
                 "Residency",
                 "Motivation"),
      labels = c(
        "Race/Ethnicity",
        "Sex",
        "Civil service exam performance",
        "Length of residency",
        "Motivation for application"
      )),
    level = word(Z, 2, sep = "_"),
    level = fct_reorder(level, estimate),
    group = "Civilian sample"
  )

g1 <- 
  gg_community_mm_pid %>%
  filter(X == "Yes") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Probability of selecting applicant",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g2 <- 
  gg_community_mm_pid %>%
  filter(X == "No") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Probability of selecting individual",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 


g3 <-
  gg_community_mm_pid %>%
  filter(X == "Difference") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Difference in selection probabilities",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g <- 
  egg::ggarrange(
    g1 + labs(x = "Democrats") + 
      theme(strip.text.x = element_text(size = 11,
                                        face = "bold"),
            panel.border = element_rect(fill = NA, colour = "black")),
    g2 + labs(x = "Non-Democrats") +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )

g

ggsave(filename = "community_conjoint_means_pid.pdf",
       g,
       width = 7.6,
       height = 3 + fixed + spacing * (nrow(g1$data)))

#### ---- Fig. S49: Differences in AMCEs by partisanship ------------------ ####

## Use this function to automate estimation process
het_amce_p <- function(x, long_data) {
  est_amce_out_x <-
    lm_robust(
      conjoint_chosen ~ sex_ + race_ + residency_ + exam_ + motivation_,
      clusters = person_id,
      data = long_data %>% filter(!!sym(x) == 1)
    ) %>%
    tidy() %>%
    filter(term != "(Intercept)") %>%
    mutate(X = "Yes")
  
  est_amce_out_nx <-
    lm_robust(
      conjoint_chosen ~ sex_ + race_ + residency_ + exam_ + motivation_,
      clusters = person_id,
      data = long_data %>% filter(!!sym(x) == 0)
    ) %>%
    tidy() %>%
    filter(term != "(Intercept)") %>%
    mutate(X = "No")
  
  est_amce_out_diff <-
    lm_robust(
      conjoint_chosen ~ !!sym(x) * sex_+!!sym(x) * race_+!!sym(x) * residency_+!!sym(x) * exam_+!!sym(x) * motivation_+!!sym(x),
      clusters = person_id,
      data = long_data
    ) %>%
    tidy() %>%
    filter(str_detect(term, ":")) %>%
    mutate(X = "Difference",
           term = str_remove(term, ".*:"))
  
  het_amce_out <-
    bind_rows(est_amce_out_x, est_amce_out_nx, est_amce_out_diff) %>%
    mutate(
      attribute = str_to_title(word(term, 1, sep = "_")),
      attribute_ref = factor(
        attribute,
        levels = c("Race",
                   "Sex",
                   "Exam",
                   "Residency",
                   "Motivation"),
        labels = c(
          "Race/Ethnicity",
          "Sex",
          "Civil service exam performance",
          "Length of residency",
          "Motivation for application"
        )
      ),
      level = word(term, 2, sep = "_"),
      level = fct_reorder(level, estimate)
    )
  return(het_amce_out)
}

gg_community_amce_pid <-
  het_amce_p(
    x = "dem_binary",
    long_data = conjoint_merged %>%
      mutate(dem_binary = as.numeric(x_pid_3_t0 == "Democrat"))
  )


g1 <- 
  gg_community_amce_pid %>%
  filter(X == "Yes") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Sub-group AMCE",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g2 <- 
  gg_community_amce_pid %>%
  filter(X == "No") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Sub-group AMCE",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 


g3 <-
  gg_community_amce_pid %>%
  filter(X == "Difference") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Difference in AMCEs",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g <- 
  egg::ggarrange(
    g1 + labs(x = "AMCE: Democrats") + 
      theme(strip.text.x = element_text(size = 10, color = "white"),
            #strip.text.x = element_blank(),
            strip.placement = "outside",
            panel.border = element_rect(fill = NA, colour = "black")),
    g2 + labs(x = "AMCE: non-Democrats") +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )

## Now convert grob into plot and add text ...
p <- ggplotify::as.ggplot(g)

pdf("community_conjoint_amces_pid.pdf",
    width = 7.6,
    height = 4 + fixed + spacing * (nrow(g1$data)))
p 
grid::grid.text("Race/ethnicity (reference: White)", 
                x = unit(0.19, "npc"), y = unit(0.985, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Sex (reference: Male)", 
                x = unit(0.19, "npc"), y = unit(0.83, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Civil service exam score (reference: Top 25%)", 
                x = unit(0.19, "npc"), y = unit(0.75, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Length of residency (reference: Does not live in city)", 
                x = unit(0.19, "npc"), y = unit(0.56, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Motivation for application (reference: Job benefits)", 
                x = unit(0.19, "npc"), y = unit(0.33, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
dev.off()

#### ---- Fig. S50: Differences in marginal means by race/ethnicity ------- ####

community_est_x0 <- list()
community_est_x1 <- list()
community_est_diff <- list()

var_list <- c(
  "sex_",
  "race_",
  "residency_",
  "exam_",
  "motivation_"
)

for (i in 1:length(var_list)) {
  col <- sym(var_list[i])
  community_est_x0[[i]] <-
    conjoint_merged %>%
    filter(x_race_cat_t0 != "White") %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i], !!col),
           X = "No")
  
  community_est_x1[[i]] <-
    conjoint_merged %>%
    filter(x_race_cat_t0 == "White") %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i], !!col),
           X = "Yes")
  
  community_est_diff[[i]] <-
    conjoint_merged %>%
    mutate(diff = as.numeric(x_race_cat_t0 == "White")) %>% 
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ diff, data = .,
      clusters = person_id
    ))) %>%
    filter(term != "(Intercept)") %>% 
    mutate(Z = paste0(var_list[i], !!col),
           X = "Difference")
  
}

gg_community_mm_race <-
  bind_rows(community_est_x0,
            community_est_x1,
            community_est_diff) %>%
  ungroup() %>%
  select(Z, X, estimate, std.error, statistic, p.value, conf.low, 
         conf.high) %>%
  mutate(
    attribute = str_to_title(word(Z, 1, sep = "_")),
    attribute_ref = factor(
      attribute,
      levels = c("Race",
                 "Sex",
                 "Exam",
                 "Residency",
                 "Motivation"),
      labels = c(
        "Race/Ethnicity",
        "Sex",
        "Civil service exam performance",
        "Length of residency",
        "Motivation for application"
      )),
    level = word(Z, 2, sep = "_"),
    level = fct_reorder(level, estimate),
    group = "Civilian sample"
  )

g1 <- 
  gg_community_mm_race %>%
  filter(X == "Yes") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Probability of selecting applicant",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g2 <- 
  gg_community_mm_race %>%
  filter(X == "No") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Probability of selecting individual",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 


g3 <-
  gg_community_mm_race %>%
  filter(X == "Difference") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Difference in selection probabilities",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g <- 
  egg::ggarrange(
    g1 + labs(x = "White respondents") + 
      theme(strip.text.x = element_text(size = 11,
                                        face = "bold"),
            panel.border = element_rect(fill = NA, colour = "black")),
    g2 + labs(x = "Non-White respondents") +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )

g

ggsave(filename = "community_conjoint_means_race.pdf",
       g,
       width = 7.6,
       height = 3 + fixed + spacing * (nrow(g1$data)))

#### ---- Fig. S51: Differences in AMCEs by race/ethnicity ---------------- #### 
gg_community_amce_race <-
  het_amce_p(
    x = "race_white",
    long_data = conjoint_merged %>%
      mutate(race_white = as.numeric(x_race_cat_t0 == "White"))
  )


g1 <- 
  gg_community_amce_race %>%
  filter(X == "Yes") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Sub-group AMCE",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g2 <- 
  gg_community_amce_race %>%
  filter(X == "No") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Sub-group AMCE",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 


g3 <-
  gg_community_amce_race %>%
  filter(X == "Difference") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Difference in AMCEs",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g <- 
  egg::ggarrange(
    g1 + labs(x = "AMCE: White respondents") + 
      theme(strip.text.x = element_text(size = 10, color = "white"),
            strip.placement = "outside",
            panel.border = element_rect(fill = NA, colour = "black")),
    g2 + labs(x = "AMCE: Non-White respondents") +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )

## Now convert grob into plot and add text ...
p <- ggplotify::as.ggplot(g)

pdf("community_conjoint_amces_race.pdf",
    width = 7.6,
    height = 4 + fixed + spacing * (nrow(g1$data)))
p 
grid::grid.text("Race/ethnicity (reference: White)", 
                x = unit(0.19, "npc"), y = unit(0.985, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Sex (reference: Male)", 
                x = unit(0.19, "npc"), y = unit(0.83, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Civil service exam score (reference: Top 25%)", 
                x = unit(0.19, "npc"), y = unit(0.75, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Length of residency (reference: Does not live in city)", 
                x = unit(0.19, "npc"), y = unit(0.56, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Motivation for application (reference: Job benefits)", 
                x = unit(0.19, "npc"), y = unit(0.33, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
dev.off()



#### ---- Fig. S52: Differences in marginal means by sex ------------------ ####

community_est_x0 <- list()
community_est_x1 <- list()
community_est_diff <- list()

var_list <- c(
  "sex_",
  "race_",
  "residency_",
  "exam_",
  "motivation_"
)

for (i in 1:length(var_list)) {
  col <- sym(var_list[i])
  community_est_x0[[i]] <-
    conjoint_merged %>%
    filter(x_female_t0 == 0) %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i], !!col),
           X = "No")
  
  community_est_x1[[i]] <-
    conjoint_merged %>%
    filter(x_female_t0 == 1) %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i], !!col),
           X = "Yes")
  
  community_est_diff[[i]] <-
    conjoint_merged %>%
    mutate(diff = as.numeric(x_female_t0 == 1)) %>% 
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ diff, data = .,
      clusters = person_id
    ))) %>%
    filter(term != "(Intercept)") %>% 
    mutate(Z = paste0(var_list[i], !!col),
           X = "Difference")
  
}

gg_community_mm_sex <-
  bind_rows(community_est_x0,
            community_est_x1,
            community_est_diff) %>%
  ungroup() %>%
  select(Z, X, estimate, std.error, statistic, p.value, conf.low, 
         conf.high) %>%
  mutate(
    attribute = str_to_title(word(Z, 1, sep = "_")),
    attribute_ref = factor(
      attribute,
      levels = c("Race",
                 "Sex",
                 "Exam",
                 "Residency",
                 "Motivation"),
      labels = c(
        "Race/Ethnicity",
        "Sex",
        "Civil service exam performance",
        "Length of residency",
        "Motivation for application"
      )),
    level = word(Z, 2, sep = "_"),
    level = fct_reorder(level, estimate),
    group = "Civilian sample"
  )

g1 <- 
  gg_community_mm_sex %>%
  filter(X == "Yes") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Probability of selecting applicant",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g2 <- 
  gg_community_mm_sex %>%
  filter(X == "No") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Probability of selecting individual",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 


g3 <-
  gg_community_mm_sex %>%
  filter(X == "Difference") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Difference in selection probabilities",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g <- 
  egg::ggarrange(
    g1 + labs(x = "Female respondents") + 
      theme(strip.text.x = element_text(size = 11,
                                        face = "bold"),
            panel.border = element_rect(fill = NA, colour = "black")),
    g2 + labs(x = "Male respondents") +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )

g

ggsave(filename = "community_conjoint_means_sex.pdf",
       g,
       width = 7.6,
       height = 3 + fixed + spacing * (nrow(g1$data)))


#### ---- Fig. S53: Differences in AMCEs by sex --------------------------- #### 
gg_community_amce_sex <-
  het_amce_p(
    x = "x_female_t0",
    long_data = conjoint_merged
  )


g1 <- 
  gg_community_amce_sex %>%
  filter(X == "Yes") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Sub-group AMCE",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g2 <- 
  gg_community_amce_sex %>%
  filter(X == "No") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Sub-group AMCE",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 


g3 <-
  gg_community_amce_sex %>%
  filter(X == "Difference") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Difference in AMCEs",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g <- 
  egg::ggarrange(
    g1 + labs(x = "AMCE: Female respondents") + 
      theme(strip.text.x = element_text(size = 10, color = "white"),
            strip.placement = "outside",
            panel.border = element_rect(fill = NA, colour = "black")),
    g2 + labs(x = "AMCE: Male respondents") +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )

## Now convert grob into plot and add text ...
p <- ggplotify::as.ggplot(g)

pdf("community_conjoint_amces_sex.pdf",
    width = 7.6,
    height = 4 + fixed + spacing * (nrow(g1$data)))
p 
grid::grid.text("Race/ethnicity (reference: White)", 
                x = unit(0.19, "npc"), y = unit(0.985, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Sex (reference: Male)", 
                x = unit(0.19, "npc"), y = unit(0.83, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Civil service exam score (reference: Top 25%)", 
                x = unit(0.19, "npc"), y = unit(0.75, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Length of residency (reference: Does not live in city)", 
                x = unit(0.19, "npc"), y = unit(0.56, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Motivation for application (reference: Job benefits)", 
                x = unit(0.19, "npc"), y = unit(0.33, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
dev.off()




#### ---- Fig. S54: Differences in marginal means by attitudes index ------ ####

## Use median split on aggregate index of trust, legitimacy, cooperation
(cut <- median(conjoint_merged$y_ltc_idx_t0, na.rm = TRUE))
conjoint_merged$ltc_binary <-
  as.numeric(conjoint_merged$y_ltc_idx_t0 >= cut)

mean(conjoint_merged$ltc_binary, na.rm = TRUE)


community_est_x0 <- list()
community_est_x1 <- list()
community_est_diff <- list()

var_list <- c(
  "sex_",
  "race_",
  "residency_",
  "exam_",
  "motivation_"
)

for (i in 1:length(var_list)) {
  col <- sym(var_list[i])
  community_est_x0[[i]] <-
    conjoint_merged %>%
    filter(ltc_binary == 0) %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i], !!col),
           X = "No")
  
  community_est_x1[[i]] <-
    conjoint_merged %>%
    filter(ltc_binary == 1) %>%
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ 1, data = .,
      clusters = person_id
    ))) %>%
    mutate(Z = paste0(var_list[i], !!col),
           X = "Yes")
  
  community_est_diff[[i]] <-
    conjoint_merged %>%
    mutate(diff = as.numeric(ltc_binary == 1)) %>% 
    group_by(!!col) %>%
    do(tidy(lm_robust(
      conjoint_chosen ~ diff, data = .,
      clusters = person_id
    ))) %>%
    filter(term != "(Intercept)") %>% 
    mutate(Z = paste0(var_list[i], !!col),
           X = "Difference")
  
}

gg_community_mm_ltc <-
  bind_rows(community_est_x0,
            community_est_x1,
            community_est_diff) %>%
  ungroup() %>%
  select(Z, X, estimate, std.error, statistic, p.value, conf.low, 
         conf.high) %>%
  mutate(
    attribute = str_to_title(word(Z, 1, sep = "_")),
    attribute_ref = factor(
      attribute,
      levels = c("Race",
                 "Sex",
                 "Exam",
                 "Residency",
                 "Motivation"),
      labels = c(
        "Race/Ethnicity",
        "Sex",
        "Civil service exam performance",
        "Length of residency",
        "Motivation for application"
      )),
    level = word(Z, 2, sep = "_"),
    level = fct_reorder(level, estimate),
    group = "Civilian sample"
  )

g1 <- 
  gg_community_mm_ltc %>%
  filter(X == "Yes") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Probability of selecting applicant",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g2 <- 
  gg_community_mm_ltc %>%
  filter(X == "No") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0.5,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Probability of selecting individual",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 


g3 <-
  gg_community_mm_ltc %>%
  filter(X == "Difference") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Difference in selection probabilities",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g <- 
  egg::ggarrange(
    g1 + labs(x = "High police legitimacy") + 
      theme(strip.text.x = element_text(size = 11,
                                        face = "bold"),
            panel.border = element_rect(fill = NA, colour = "black")),
    g2 + labs(x = "Low police legitimacy") +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 11),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )

g

ggsave(filename = "community_conjoint_means_ltc.pdf",
       g,
       width = 7.6,
       height = 3 + fixed + spacing * (nrow(g1$data)))


#### ---- Fig. S55: Differences in AMCEs by attitudes index --------------- ####
gg_community_amce_ltc <-
  het_amce_p(
    x = "ltc_binary",
    long_data = conjoint_merged
  )


g1 <- 
  gg_community_amce_ltc %>%
  filter(X == "Yes") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 22,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Sub-group AMCE",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    strip.text.x = element_text(hjust = 0),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )

g2 <- 
  gg_community_amce_ltc %>%
  filter(X == "No") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 21,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Sub-group AMCE",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  ) 


g3 <-
  gg_community_amce_ltc %>%
  filter(X == "Difference") %>%
  ggplot(.,
         aes(x = estimate,
             y = level)) +
  geom_vline(
    xintercept = 0,
    color = "black",
    lty = 2,
    size = 0.2
  ) +
  geom_errorbarh(
    aes(xmin = conf.low, xmax = conf.high),
    height = 0,
    color = "black",
    position = position_dodgev(height = .5),
    size = 1
  ) +
  geom_point(
    size = 2.5,
    pch = 24,
    color = "white",
    fill = "black",
    position = position_dodgev(height = .5)
  ) +
  facet_col(attribute_ref ~ ., scales = "free_y", space = "free") +
  labs(x = "Difference in AMCEs",
       y = "",
       caption = "") +
  theme_tufte() +
  theme(
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(size = 0.25),
    strip.text = element_text(size = 11),
    axis.line.x = element_line(size = 0.25),
    axis.title.x = element_text(size = 10),
    axis.text.y = element_text(size = 11),
    plot.margin = unit(c(
      t = 0.1,
      r = 0.1,
      b = -0.5,
      l = 0.1
    ), "lines"),
    legend.position = "bottom",
    legend.justification = "left",
    legend.margin = margin(0, 0, 0, 0),
    legend.box.margin = margin(-10, -10, -10, -10)
  )


g <- 
  egg::ggarrange(
    g1 + labs(x = "AMCE: High police legitimacy") + 
      theme(strip.text.x = element_text(size = 10, color = "white"),
            strip.placement = "outside",
            panel.border = element_rect(fill = NA, colour = "black")),
    g2 + labs(x = "AMCE: Low police legitimacy") +
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    g3 + labs(x = "Difference")  + 
      theme(
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        strip.text.x = element_text(color = "white", size = 10),
        panel.border = element_rect(fill = NA, colour = "black")
      ),
    nrow = 1
  )

## Now convert grob into plot and add text ...
p <- ggplotify::as.ggplot(g)

pdf("community_conjoint_amces_ltc.pdf",
    width = 7.6,
    height = 4 + fixed + spacing * (nrow(g1$data)))
p 
grid::grid.text("Race/ethnicity (reference: White)", 
                x = unit(0.19, "npc"), y = unit(0.985, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Sex (reference: Male)", 
                x = unit(0.19, "npc"), y = unit(0.83, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Civil service exam score (reference: Top 25%)", 
                x = unit(0.19, "npc"), y = unit(0.75, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Length of residency (reference: Does not live in city)", 
                x = unit(0.19, "npc"), y = unit(0.56, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
grid::grid.text("Motivation for application (reference: Job benefits)", 
                x = unit(0.19, "npc"), y = unit(0.33, "npc"), 
                just = "left",
                gp = gpar(fontsize = 11, col = "black", fontfamily = "serif",
                          fontface = "bold"))
dev.off()



